#include <iostream>
#include <stdlib.h>
#include <time.h>
#include <math.h>

#define PI 3.14159


using namespace std;

int slice_num = 0;
double t = 0;
double a = 0;

void initialize_rand(){
	// Here we use time as seed
	srand((unsigned long)time(NULL));
}

//We need to know that this rand generator is not very good 
double normalization_rand(){
	return 1.0 * rand() / RAND_MAX;
}

double get_f(double *x){
	return x[0] * x[0];
}


double* Initialize(){
	double *x = new double[slice_num];
	for (int i = 0; i < slice_num; i++){
		x[i] = 0.0;
	}
	return x;
}

double get_delta(double n){
	return (normalization_rand() - 0.5) * n;
}

double get_potential(double x){
	return 0.5 * x * x;
}

int main(int argc, char *argv[]){
	
	t = atof(argv[1]);
	slice_num = 20;
	a = (1.0 / t ) / slice_num;
	//cout << t << endl;
	//cout << a << endl;
	long n = 4000000;
	long pre_n = 20000;
	long accept_test = 1000;
	long accept = 0;
	double energy = 0;
	double delta_E;
	double x_trial;
	double delta;
	double ratio;
	double *x = Initialize();
	initialize_rand();
	
	double coe = 1.0;
	while (1){
		accept = 0;
    	for (long i = 0; i < accept_test; i++){
        	for (int j = 0; j < slice_num; j++){
            	delta = get_delta(coe);
            	x_trial = x[j] + delta;
            	if ( j == 0){
                	delta_E = 1.0 / a * delta * (x_trial + x[j] - x[j + 1] - x[slice_num - 1]) + a * (get_potential(x_trial) - get_potential(x[j]));
            	}
            	else if ( j == (slice_num - 1)){
                	delta_E = 1.0 / a * delta * (x_trial + x[j] - x[0] - x[j - 1]) + a * (get_potential(x_trial) - get_potential(x[j]));
            	}
            	else{
                	delta_E = (1.0 / a) * delta * (x_trial + x[j] - x[j + 1] - x[j - 1]) + a * (get_potential(x_trial) - get_potential(x[j]));
            	}	
            	if (delta_E <= 0 || normalization_rand() <= exp(-1.0 * delta_E)){
                	x[j] = x_trial;
					accept++;
            	}
        	}
    	}
		ratio = 1.0 * accept / (accept_test * slice_num);
		if (ratio > 0.48){
			coe = 1.05 * coe;
		}
		else if (ratio < 0.41){
			coe = 0.95 * coe;
		}
		else
			break;
	}

	//cout << ratio << endl;

	for (long i = 0; i < pre_n; i++){
		for (int j = 0; j < slice_num; j++){
			delta = get_delta(coe);
            x_trial = x[j] + delta;
            if ( j == 0){
                delta_E = 1.0 / a * delta * (x_trial + x[j] - x[j + 1] - x[slice_num - 1]) + a * (get_potential(x_trial) - get_potential(x[j]));
            }
            else if ( j == (slice_num - 1)){
                delta_E = 1.0 / a * delta * (x_trial + x[j] - x[0] - x[j - 1]) + a * (get_potential(x_trial) - get_potential(x[j]));
            }
            else{
                delta_E = (1.0 / a) * delta * (x_trial + x[j] - x[j + 1] - x[j - 1]) + a * (get_potential(x_trial) - get_potential(x[j]));
			}
			if (delta_E <= 0 || normalization_rand() <= exp(-1.0 * delta_E)){
                x[j] = x_trial;
			}
		}
	}
	
	accept = 0;
	energy = 0;
	double temp_value;
	for (long i = 0; i < n; i++){
		double temp = 0;
		for (int j = 0; j < slice_num; j++){
			delta = get_delta(coe);
			x_trial = x[j] + delta;
			if (j == 0){
				delta_E = 1.0 / a * delta * (x_trial + x[j] - x[j + 1] - x[slice_num - 1]) + a * (get_potential(x_trial) - get_potential(x[j]));			
			}
			else if (j == (slice_num - 1)){
				delta_E = 1.0 / a * delta * (x_trial + x[j] - x[0] - x[j - 1]) + a * (get_potential(x_trial) - get_potential(x[j]));
			}
			else{
				delta_E = (1.0 / a) * delta * (x_trial + x[j] - x[j + 1] - x[j - 1]) + a * (get_potential(x_trial) - get_potential(x[j]));
				//delta_E = a * ( 0.5 * ((x[j + 1] - x_trial) * (x[j + 1] - x_trial) + (x_trial - x[j - 1]) * (x_trial - x[j - 1])) / (a * a) + get_potential(x_trial)) - a * ( 0.5 * ((x[j + 1] - x[j]) * (x[j + 1] - x[j]) + (x[j] - x[j - 1]) * (x[j] - x[j - 1])) / (a * a) + get_potential(x[j]));
			}

			//cout << delta_E << endl;;	
			if (delta_E <= 0 || normalization_rand() <= exp(-1.0 * delta_E)){
				x[j] = x_trial;
				//cout << delta << endl;
				accept++;
			}
			
			temp += (x[j] * x[j]);
		}
		temp_value += temp;
		if ((i + 1 )  % pre_n == 0){
			energy += temp_value;
			temp_value = 0;
			cout << t << "\t" << 1.0 * accept / (slice_num * (i + 1)) << "\t" << energy / ((i + 1.0) * (slice_num))  << endl;
		}
	}
	//cout << t << "\t" << 1.0 * accept / (slice_num * n) << "\t" << energy / ((1.0 * n) * (slice_num)) << endl;
	/*
	for (int j = 0; j < slice_num; j++){
		cout << j << "\t" << x[j] << endl;
	}
	*/
	free(x);
	return 0;
}
